;+
; NAME:
;       BGPS_FFORE
;
; PURPOSE:
;       Model (using chi-square minimization) the foreground fraction
;       of 8-micron light along a line of sight to an IRDC
;
; CATEGORY:
;       IRDC modeling
;
; CALLING SEQUENCE:
;       BGPS_FFORE [,Td][,MODELNO=modelno]
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       Td      -- Dust temperature [Default: fit an ensemble
;                  temperature to the entire data set]
;       CASENO  -- Case Number of the case(s) to be run. May be
;                  scalar or vector.  Cases are as follows:
;                  1) stars + gas, no veil
;                  2) stars + veil, no gas
;                  3) gas + veil, no stars
;                  4) veil + gas + stars
;       FN      -- Name of the RB3 structure IDL save file to use
;                  [Default: bgps_rb3.sav]
;
; KEYWORD PARAMETERS:
;       CLEM -- Use the Clemens (1985) rotation curve instead of the
;               Reid et al (2009) curve.
;       IAU  -- Use the IAU Standard rotation curve instead of the
;               Reid et al (2009) curve.
;       If alternate rotation curve is used, a suffix in the output
;       filename will be appended.
;
; OUTPUTS:
;       Writes ISL sace file case?_irdc_Td???.sav
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       FFORE_BLOCK
;
; MODIFICATION HISTORY:
;
;       Created:  03/07/11, TPEB -- Initial codified version.
;       Modified: 03/08/11, TPEB -- Name change to model_ffore & added
;                                   option to specify dust
;                                   temperature.
;       Modified: 03/28/11, TPEB -- Added models and changed array and
;                                   model indexing to be more
;                                   intuitive.  Also, added
;                                   functionality for returning WHICH
;                                   sources are at the far kdist.
;                                   ALSO: fixed bug in calculating
;                                   sig_f -- was set for only Td=15K.
;       Modified: 04/05/11, TPEB -- Adjustments for grid spacing and
;                                   ranges.
;       Modified: 04/29/11, TPEB -- Added driver for Model 13.
;       Modified: 05/02/11, TPEB -- Changed model selection to be a
;                                   CASE statement, so don't
;                                   have to keep commenting out stuff.
;       Modified: 05/17/11, TPEB -- Added fuctionality to fit for Td
;                                   if input Td not specified.  (Plus
;                                   documentation and code cleanup.)
;       Modified: 05/23/11, TPEB -- Moved calculation of (R,Z) vectors
;                                   to the new routine omni_lbd2rz.pro,
;                                   which takes into account
;                                   inaccuracies in the definition of
;                                   Galactic coordinates.
;       Modified: 05/24/11, TPEB -- Placed calculation of (R,Z)
;                                   vectors and kinematic distances
;                                   into the main routine and passed
;                                   to the LOOPIE via a structure in
;                                   the FFORE_BLOCK.  This should
;                                   greatly speed calculations since
;                                   previously these were each
;                                   calculated for every source for
;                                   every value of (A,B).
;       Modified: 06/29/11, TPEB -- Added functionality for choosing
;                                   rotation curves other than the
;                                   Reid et al. (2009) model assumed
;                                   up to this point.  Added the IAU
;                                   standard curve and the Clemens
;                                   (1985) curve.  The choice of
;                                   rotation curve affects the
;                                   kinematic distances calculated for
;                                   each object.
;       Modified: 07/22/11, TPEB -- Added option to run on the BGPS
;                                   data set contained in the IDL save
;                                   file irdc_dist_model/bgsp_rb3.sav
;       Created:  08/10/11, TPEB -- Cut a branch into the routine
;                                   BGPS_FFORE.pro, a routine
;                                   specifically for modeling the
;                                   foreground fraction using BGPS
;                                   data.
;       Modified: 08/12/11, TPEB -- Moved ALL repetitive calculations
;                                   to the driver routine, so be
;                                   shared with the other routines via
;                                   the FFORE_BLOCK.
;       Modified: 08/17/11, TPEB -- Added option for using different
;                                   RB3 save file (for training sets
;                                   or other reasons).
;       Modified: 08/25/11, TPEB -- Changed output method from an
;                                   ASCII text file containing each
;                                   point along the trajectory to an
;                                   IDL save file containing the
;                                   best-fit values along with the
;                                   various statistics useful for
;                                   later analysis.
;       Modified: 08/31/11, TPEB -- Moved output IDL save files to
;                                   subdirectory irdc_dist_model/data/
;       Modified: 09/06/11, TPEB -- Ensured all BGPS_FFORE_MODEL
;                                   calculations are being done with
;                                   double().
;       Modified: 09/06/11, TPEB -- Added subroutine STAR_DENSITY to
;                                   place that calculation only in one
;                                   place (rather than copies
;                                   everywhere). 
;       Modified: 12/12/11, TPEB -- Adding limits on the size of step
;                                   allowed by MPFIT for the
;                                   parameters in the fit.  This is to
;                                   (partially?) compensate for the
;                                   (uncontrollable) effect of placing
;                                   sources at N/F based on chisq
;                                   values.
;       Modified: 01/09/12, TPEB -- Added NEARFAR to the IDL save
;                                   file.
;       Modified: 03/09/12, TPEB -- Testing a modification that would
;                                   force f_data to be >= 0. (BAD
;                                   IDEA!)
;       Modified: 03/30/12, TPEB -- Small modification to MPFIT
;                                   parinfo definitions to make it
;                                   easier for a human to read.
;       Modified: 05/02/12, TPEB -- Fixed normalization in
;                                   STAR_DENSITY to accurately return
;                                   Sigma(R_sun) = 30 M_sun / pc^2.
;       Modified: 09/04/12, TPEB -- Adjusted parametric models to have
;                                   a Gaussian vertical distribution
;                                   rather than exponential.
;
;-

;;=================================================================
;; Stellar disk density function
;;=================================================================
FUNCTION STAR_DENSITY, R, Z
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  
  R0 = double(!MW.R0)
  R  = double(R)
  Z  = double(Z)
  
  ;; The 800 = 2/3 * (H1 + H2) is the normalization
  ;; sigma = 30 M_sun / pc^2 at the solar circle
  
  rho = 30.d * (exp( -(R-R0)/2600.d - abs(Z)/300.d) + $
                0.12d * exp( -(R-R0)/3600.d - abs(Z)/900.d)) / (800.d)
  
  RETURN, rho
END

;;=================================================================
;; Single, generic model w/ 3 parameters
;;=================================================================
FUNCTION BGPS_FFORE_MODEL, dist, A, B, C, NG, rho
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  rho = NG * (B * rho_hi + rho_h2) + (C * rho_star)
  
  ;; Integration limits
  ilt = WHERE(d LE dist, nlt)
  IF nlt GT 0 THEN $
     RETURN,A + (TOTAL(rho[ilt]) / TOTAL(rho)) ELSE $
        RETURN, A
  
END


;;=================================================================
;; Single, generic model w/ 3 parameters, alternate version
;;=================================================================
FUNCTION BGPS_FFORE_MODEL2, dist, A, B, C, rho
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  rho = (B * rho_hi) + (C * rho_h2) + rho_star
  
  ;; Integration limits
  ilt = WHERE(d LE dist, nlt)
  IF nlt GT 0 THEN $
     RETURN,A + (TOTAL(rho[ilt]) / TOTAL(rho)) ELSE $
        RETURN, A
  
END


;;=================================================================
;; Simple, parametric model, used to test constraints on rho_gas
;;=================================================================
FUNCTION BGPS_FFORE_SIMPLE, dist, A, L, H, rho
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  
  R0 = double(!MW.R0)
  R  = double(R)
  Z  = double(Z)
  
  rho = exp( -(R-R0)/L)* exp(-Z*Z/H/H/2.d)
  
  ;; Integration limits
  ilt = WHERE(d LE dist, nlt)
  IF nlt GT 0 THEN $
     RETURN,A + (TOTAL(rho[ilt]) / TOTAL(rho)) ELSE $
        RETURN, A
  
END


;;=================================================================
;; Simple, parametric model, used to test constraints on rho_gas, PLUS
;;   emission to simulate Molecular Ring
;;=================================================================
FUNCTION BGPS_FFORE_SIMPLE2, dist, A, L, H, GA, GW, rho
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  
  R0 = double(!MW.R0)
  R  = double(R)
  Z  = double(Z)
  
  rho = (exp( -(R-R0)/L) + GA *exp(-(R-4500.d)^2 / (2.*GW*GW)  )) * $
        exp(-Z*Z/H/H/2.d) * ( exp( (R-3600.d)/500.d ) < 1.d )
  
  ;; Integration limits
  ilt = WHERE(d LE dist, nlt)
  IF nlt GT 0 THEN $
     RETURN,A + (TOTAL(rho[ilt]) / TOTAL(rho)) ELSE $
        RETURN, A
  
END


;;=================================================================
;; Function for use with TNMIN to return the chisq value based on
;; an input temperature
;;=================================================================
FUNCTION BGPS_FFORE_FITT, p, KD=kd, A=A, B=B, C=C, NG=NG, RHO=rho,$
                          OBJNO=i, RESID=resid, _EXTRA=extra
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  TEMP_DEPENDENT_VARIABLES, p[0], OBJNO=i
  f_model = BGPS_FFORE_MODEL(kd, A, B, C, NG)
  
  resid = (f_model - f_data[i]) / sig_f[i]
  
  RETURN,resid * resid
END


;;=================================================================
;; Model - I_MIR scaling
;;=================================================================
FUNCTION BGPS_FFORE_IMIR_SCALE, P, IMIR=imir, MODEL=model
  
  ind = WHERE(imir NE 0, nind)
  resid = imir - (model * P[0])
  RETURN, resid[ind]
END

;;=================================================================
;;=================================================================
;; LOOPIE FUNCTION FOR USE WITH MPFIT -- a way of bypassing the
;;   lengthy gridding operation
;;=================================================================
;; Function that does the looping over sources for a given A & B
FUNCTION BGPS_FFORE_LOOPIE, P, VERBOSE=verbose,NEARFAR=nearfar, $
                            FFORE_VEC=ffore_vec, FDIFF=fdiff, $
                            CONSTRAIN=constrain
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  P = double(P)
  A  = P[0]
  B  = P[1]
  C  = P[2]
  NG = P[3]
  Td = P[4]
  mt = P[5]
  GW = P[6]
  
  resid     = dblarr(n)
  ffore_vec = dblarr(n)
  fdiff     = dblarr(n)
  nearfar   = bytarr(n)
  constrain = bytarr(n)
  
  ;; Choice here.  If NOT fitting temperature, calculate model
  ;;    ffore and chi^2_i as usual.  Else...
  IF do_Tfit THEN TEMP_DEPENDENT_VARIABLES, Td
  
  chisq = 0.d
  FOR i=0L, n-1 DO BEGIN
     
     ;; Place LPSTR elements into regular variables for use in this LOOPIE
     R        = lpstr[i].R
     Z        = lpstr[i].Z
     kdn      = lpstr[i].kdn
     kdf      = lpstr[i].kdf
     rho_hi   = lpstr[i].rho_hi
     rho_h2   = lpstr[i].rho_h2
     rho_star = lpstr[i].rho_star
     
     ;; Select where to get model
     CASE mt OF
        0: BEGIN
           f_model_n = BGPS_FFORE_MODEL(kdn, A, B, C, NG)
           f_model_f = BGPS_FFORE_MODEL(kdf, A, B, C, NG)
        ENDCASE
        1: BEGIN
           f_model_n = BGPS_FFORE_SIMPLE(kdn, A, B, C)
           f_model_f = BGPS_FFORE_SIMPLE(kdf, A, B, C)
        ENDCASE
        2: BEGIN
           f_model_n = BGPS_FFORE_MODEL2(kdn, A, B, C)
           f_model_f = BGPS_FFORE_MODEL2(kdf, A, B, C)
        ENDCASE
        3: BEGIN
           f_model_n = BGPS_FFORE_SIMPLE2(kdn, A, B, C, NG, GW)
           f_model_f = BGPS_FFORE_SIMPLE2(kdf, A, B, C, NG, GW)
        ENDCASE
        ELSE:  message,'Model Type '+string(mt,format="(I0)")+' not valid.'
     ENDCASE
     
     dcn = (f_model_n - f_data[i])^2 / sig_f[i]^2
     dcf = (f_model_f - f_data[i])^2 / sig_f[i]^2
     residn = (f_model_n - f_data[i]) / sig_f[i]
     residf = (f_model_f - f_data[i]) / sig_f[i]
     
     ;; Marginalize over NEAR / FAR
     delch  = dcn < dcf
     fdiff[i] = abs(dcn - dcf) / delch
     ;; Well-constrained if difference between chisq is large, if the
     ;; object is at d_tan (fdiff[i] = 0), or if it is within 0.5 kpc
     ;; of the tangent distance.
     constrain[i] = ((fdiff[i] GE 1.0) OR (fdiff[i] EQ 0)) OR $
                    (abs(kdf-kdn) LE 1.d3)
     nearfar[i]  = dcf LE dcn
     chisq += delch
     resid[i]     = nearfar[i] ? residf : residn
     ffore_vec[i] = nearfar[i] ? f_model_f : f_model_n
     
  ENDFOR
  
  RETURN,resid
  
END


;; =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
;; MAIN Driver routine =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
;; =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=

PRO BGPS_FFORE, Td, CASENO=caseno, CLEM=clem, IAU=iau, FN=fn, DOIMIR=doimir
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  FORWARD_FUNCTION IMIR_MOSAIC_MODEL
  
  ;; Get input for which model(s) to run.
  IF n_elements( caseno ) EQ 0 THEN BEGIN
     message,'Enter case number to run (1-4):',/inf
     modstr = ''
     read,modstr
     caseno = fix(modstr)
  ENDIF
  donmod = n_elements(caseno)
  do_imir_junk = KEYWORD_SET( doimir )
  
  ;; Choose rotation curve to use
  IF KEYWORD_SET( clem ) THEN rotc = 2b ELSE $
     rotc = KEYWORD_SET( iau ) ? 1b : 0b
  
  ;;=====================================================================
  ;; Print system time for timing purposes
  print,'Start Time:' & spawn,'date'
  
  ;; Read in the IDL save file RB3
  IF n_elements(fn) EQ 0 THEN fn = 'bgps_rb3.sav'
  IDIR = './irdc_dist_model/'
  restore,IDIR+fn,/ver
  subset = strmid(fn,strpos(fn,'_rb3')+4,-strpos(fn,'_rb3')-4+strpos(fn,'.sav'))
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  corr = R0 / 8.5d3
  
  ;; Do some work, and place this info into the structure!
  n       = n_elements(rb3)
  tau     = dblarr(n)
  f_data  = dblarr(n)
  sig_f   = dblarr(n)
  
  ;; Do the re-usable calculations
  GENERATE_LPSTR, ROTC=rotc, SUFF=suff
  
  ;; Check if doing fitting over temperature... take appropriate action.
  do_Tfit = (n_elements(Td) EQ 0) ? 1b : 0b
  IF ~ do_Tfit THEN TEMP_DEPENDENT_VARIABLES, Td
  
  message,'Starting the fitting...',/inf
  
  ;;================*****===========
  ;; Enter present number of models (or many more) -- only generates
  ;;   an array of filenames.
  
  nmods = 7
  dfn = strarr(nmods+1)
  infix = do_Tfit ? '_irdc_TdFIT' : '_irdc_Td'+string(Td,format="(I0,'K')")
  
  FOR j=1, nmods DO $
     dfn[j] = IDIR+'data/case'+string(j,format="(I0)")+infix+subset+suff+'.sav'
  
  
  ;;======================================================
  ;; CASE Statement for choosing cases
  
  FOR k = 0L, donmod-1 DO BEGIN
     
     ;; Set up the general framework for the MPFIT scheme here.
     
     parinfo = replicate( $
               create_struct('VALUE',0.d,'FIXED',0b,$
                             'LIMITED',[1b,1b],'LIMITS',[0.d,1.d],$
                             'PARNAME','','MPPRINT',1b,'MPMAXSTEP',0.d), 7)
     parinfo.parname = ['Veil Amount','Relative HI','Stellar SDn',$
                        'Neutral Gas','Dust Temp  ','Model Type ','']
     parinfo[1].limits  = [0.d,100.d]
     parinfo[2].limits  = [0.d,100.d]
     parinfo[4].limits  = [10.d,40.d]
     parinfo[5].limited = [0b,0b]
     parinfo[5].fixed = 1b
     parinfo[6].fixed = 1b
     parinfo[6].mpprint = 0b
     
     ;; Move Td by no more than 1K per step, others, too.
     parinfo.mpmaxstep = [0.02d,5.d,5.d,0.d,1.d,0.d,0.d]
     
     IF ~ do_Tfit THEN BEGIN
        parinfo[4].fixed = 1b
        parinfo[4].value = Td
     ENDIF ELSE BEGIN
        parinfo[4].value = 20.d
     ENDELSE
     
     
     CASE caseno[k] OF
        
        ;; First Case ====================================================
        10: BEGIN
           message,'Case I:  (stars + gas, no veil)',/inf
           savfn = dfn[1]
           
           ;; Set parinfo for this case
           parinfo[0:3].value   = [0.d, 2.5d, 1.d, 1.d]
           parinfo[0:3].fixed   = [1b,0b,0b,1b]
           parinfo[0:3].mpprint = [0b,1b,1b,1b]
           
        ENDCASE
        
        ;; Second Case ====================================================
        20: BEGIN
           message,'Case II:  (stars + veil, no gas)',/inf
           savfn = dfn[2]
           
           ;; Set parinfo for this case
           parinfo[0:3].value   = [0.10d, 0.d, 1.d, 0.d]
           parinfo[0:3].fixed   = [0b,1b,0b,1b]
           parinfo[0:3].mpprint = [1b,0b,1b,1b]
           parinfo[2].limits = [0.d,1.d]   ;; With no gas, no need to go to 100
           
        ENDCASE
        
        ;; Third Case ====================================================
        30: BEGIN
           message,'Case III:  (gas + veil, no stars)',/inf
           savfn = dfn[3]
           
           ;; Set parinfo for this case
           parinfo[0:3].value   = [0.10d, 2.5d, 0.d, 1.d]
           parinfo[0:3].fixed   = [0b,0b,1b,1b]
           parinfo[0:3].mpprint = [1b,1b,0b,1b]
           
        ENDCASE
        
        ;; Fourth Case ====================================================
        40: BEGIN
           message,'Case IV:  (veil + gas + stars)',/inf
           savfn = dfn[4]
           
           ;; Set parinfo for this case
           parinfo[0:3].value   = [0.10d, 2.5d, 1.d, 1.d]
           parinfo[0:3].fixed   = [0b,0b,0b,1b]
           parinfo[0:3].mpprint = [1b,1b,1b,1b]
           
        ENDCASE
        
        ;; Latecomer, 4th case with rho_star fixed
        3: BEGIN
           message,'Case V:  (veil + gas + stars), rho_star fixed',/inf
           savfn = dfn[3]
           
           ;; Reset just about everything in parinfo
           parinfo[0:3].parname = ['Veil Amount','Relative HI','Relative H2',$
                                   'Neutral Gas']
           
           ;; Set parinfo for this case
           parinfo[0:3].value   = [0.10d, 2.5d, 1.d, 1.d]
           parinfo[0:3].fixed   = [0b,0b,0b,1b]
           parinfo[0:3].mpprint = [1b,1b,1b,1b]
           parinfo[5].value = 2.
           
        ENDCASE
        
        ;; Latecomer: Simplified parametric case
        1: BEGIN
           message,'Case VI:  (parametric exponential disk)',/inf
           savfn = dfn[1]
           
           ;; Reset just about everything in parinfo
           parinfo.parname = ['Veil Amount','Rad. Scale ','Vert. Scale',$
                              '           ','Dust Temp  ','Model Type ',$
                              '           ']
           
           ;; Set parinfo for this case
           parinfo[1].limits  = [500.d,10000.d]
           parinfo[2].limits  = [10.d,1500.d]
           parinfo[3].limits  = [0.d,500.d]
           parinfo[0:3].value   = [0.0d, 2500.d, 100.d, 10.d]
           parinfo[0:3].fixed   = [0b,0b,0b,1b]
           parinfo[0:3].mpprint = [1b,1b,1b,0b]
           parinfo[5].value = 1.

           ;; Move Td by no more than 1K per step, others, too.
           parinfo.mpmaxstep = [0.01d,200.d,50.d,0.d,1.d,0.d,0.d]
           
        ENDCASE

        ;; Latecomer: Simplified parametric case + goodies
        2: BEGIN
           message,'Case VII:  (parametric exponential disk + goodies)',/inf
           savfn = dfn[2]
           
           ;; Reset just about everything in parinfo
           parinfo.parname = ['Veil Amount','Rad. Scale ','Vert. Scale',$
                              'Gauss Amp  ','Dust Temp  ','Model Type ',$
                              'Gauss Width']
           
           ;; Set parinfo for this case
           parinfo[1].limits  = [500.d,10000.d]
           parinfo[2].limits  = [10.d,1500.d]
           parinfo[3].limits  = [0.d,500.d]
           parinfo[0:3].value   = [0.0d, 2500.d, 100.d, 3.d]
           parinfo[0:3].fixed   = [0b,0b,0b,0b]
           parinfo[0:3].mpprint = [1b,1b,1b,1b]
           parinfo[5].value = 3.
           parinfo[6].fixed = 0b
           parinfo[6].value = 250.d
           parinfo[6].limits = [100.,2000.]
           parinfo[6].mpprint = 1b

           ;; Move Td by no more than 1K per step, others, too.
           parinfo.mpmaxstep = [0.01d,200.d,50.d,1.d,1.d,0.d,100.d]
           
        ENDCASE
        
        ELSE: message,'Improper case number specified...'
     ENDCASE
     
     functargs = {VERBOSE:0b}
     res = MPFIT('BGPS_FFORE_LOOPIE',parinfo.value, $
                 ERRMSG=errmsg, PARINFO=parinfo, STATUS=status, $
                 FUNCTARGS=functargs, COVAR=covar, BESTNORM=bestnorm, $
                 DOF=dof, NPEGGED=npegged)
     
     ;;========================================================================
     ;;========================================================================
     ;; Crazy-ass IMIR fitting section
     ;;
     IF do_imir_junk THEN BEGIN
        ;; At this point, run IMIR_MOSAIC
        tempfn = IDIR+STRMID(fn, 0, STRPOS(fn,'.sav'))+'_bgpsffore.sav'
        save,rb3,filename=tempfn,/ver
        print,tempfn
        IMIR_MOSAIC, SUBSET='_bgpsffore', SUCCESS=success
        IF ~ success THEN message,'Sadnesses...'
        imirm = readfits('./local/output/Imir_mosaic'+subset+'_bgpsffore.fits',$
                         mhd)
        
        ;; Re-pixelate to 10x coarser image for this step
        smim = fltarr(fix(45./(72./3600.))+2, fix(2.5/(72./3600.))+2)
        mkhdr,refhd,smim
        smsz = long(size(smim,/DIM))
        make_astr,astr,CD=[[-0.0199999986216d,0.],[0.,0.0199999986216d]],$
                  CTYPE=['GLON-CAR','GLAT-CAR'],CRPIX=smsz/2L,CRVAL=[42.5d,0.d]
        putast,refhd,astr
        hastrom,imirm,mhd,refhd
        
        modimir = IMIR_MOSAIC_MODEL(imirm,res,parinfo[5].value,astr,mhd)
        modimir *= (MPFIT('BGPS_FFORE_IMIR_SCALE', [1.], $
                          PARINFO={value:1.d,limited:[1,0],limits:[0.,0.], $
                                   parname:'Scale Factor'},$
                          functargs={IMIR:imirm, MODEL:modimir},$
                          ERRMSG=errmsg,STATUS=status))[0]
        writefits,'./local/output/Imir_mosaic'+subset+'_bgpsffore_model.fits',$
                  modimir,mhd
        
        ;; Plotting
        undefine,range
        myps,'./irdc_dist_model/analysis_plots/imir_vs_model.eps',xsize=15,$
             ct=13,/cmyk
        multiplot,[1,2],ygap=0.04,/doxaxis,/doyaxis

        plotmap,imirm,mhd,RANGE=range
        loadct,13,/silent
        cgColorbar,RANGE=range,title='I!dMIR!n Intensity [MJy sr!u-1!n]',$
                   charsize=1.0,position=[0.10,0.87,0.90,0.92]
        multiplot,/doyaxis,/doxaxis
        plotmap,modimir,mhd,RANGE=range
        
        cgText,/norm,0.85,0.6,'Milky Way',charsize=1.2,color='Blue Violet'
        cgText,/norm,0.85,0.22,'Model, Case 2',charsize=1.2,color='Blue Violet'
        
        myps,/done,/mp
        undefine,imirm
     ENDIF
     ;;
     ;; End of crazy-ass IMIR fitting section
     ;;========================================================================
     ;;========================================================================
     
     ;; Just run the loopie on the result to get nearfar
     junk = BGPS_FFORE_LOOPIE( res, NEARFAR=nearfar, FDIFF=fdiff, $
                               CONSTRAIN=constrain, FFORE_VEC=ffore_vec)
     
     nfar = fix(total(nearfar))
     print,'NFAR: ',nfar
     correl = covarnorm(covar)
     print,status,errmsg
     print,'Best-Fit values are: ',res
     print,covar
     print,'CHISQ: ',bestnorm,bestnorm/double(dof+npegged)
     
     ;; More stuff to save...
     kdn = lpstr.kdn
     kdf = lpstr.kdf
     fitdist = nearfar*kdf + (~nearfar)*kdn
     
     save,res,bestnorm,dof,npegged,covar,correl,nfar,nearfar,ffore_vec,fdiff,$
          constrain,kdn,kdf,fitdist,filename=savfn,/ver
     
     plusm = string("261B);"
     
     print,string(res[4],plusm,sqrt(covar[4,4]),format=$
                  "('Temperature: ',F0.1,' ',A0,' ',F0.1)")
     
     print,correl
     
  ENDFOR
  
  ;;======================================
  ;; End of Code
  eoc:
  print,'End Time:' & spawn,'date'
  
END
